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Abstract 

Safe  and  reliable  operation  of  a  fuel  cell  requires  proper  management  of  the  water  and  heat  that  are  produced  as  by-products.  Most  of 
the  current  models  for  the  cell  used  for  an  analysis  of  the  fuel  cell  system  are  based  on  the  empirical  polarization  curve  and  neglect  the 
dynamic  effects  of  water  concentration,  temperature  and  reactant  distribution  on  the  characteristics.  The  new  model  proposed  in  this  paper 
is  constructed  upon  the  layers  of  a  cell,  taking  into  account  the  following  factors:  (1)  dynamics  in  temperature  gradient  across  the  fuel  cell; 
(2)  dynamics  in  water  concentration  redistribution  in  the  membrane;  (3)  dynamics  in  proton  concentration  in  the  cathode  catalyst  layer;  (4) 
dynamics  in  reactant  concentration  redistribution  in  the  cathode  GDL.  Simulations  have  been  performed  to  analyze  the  effects  of  load  currents 
on  the  behaviors  of  the  fuel  cell.  In  the  future,  the  fuel  cell  model  will  be  extended  to  a  stack  model  and  integrated  with  system  models.  All 
of  the  models  will  be  implemented  on  a  real  time  system  that  optimizes  the  computation  time  by  a  parallelization  of  solvers,  which  provides 
an  environment  to  analyze  the  performance  and  optimize  design  parameters  of  the  PEM  fuel  cell  system  and  components. 

©  2005  Published  by  Elsevier  B.Y. 
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1.  Introduction 

The  PEM  fuel  cell  is  a  strong  candidate  for  use  as  an 
alternative  power  source  in  future  vehicle  and  power  condi¬ 
tioning  applications.  The  effects  of  electric  loads  on  tem¬ 
perature,  water  in  the  stack  and  reactants  are  crucial  is¬ 
sues  that  must  be  considered  for  the  optimum  design  of 
fuel  cell  powered  systems.  Currently,  fuel  cell  stack  mod¬ 
els  are  being  employed  to  analyze  these  effects.  However, 
the  simulation  results  do  not  incorporate  either  the  dynamic 
or  transient  aspects  of  the  fuel  cell  system  in  operating 
environments. 

As  a  matter  of  fact,  the  dynamic  power  output  and  ef¬ 
ficiency  profile  of  a  PEMFC  is  strongly  influenced  by  the 
variation  of  the  temperature,  reactant  and  product  transfer  in 
the  fuel  cell  caused  by  a  current  load. 
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Firstly,  the  temperature  significantly  affects  the  perfor¬ 
mance  of  a  fuel  cell  by  influencing  the  water  removal  and 
reactants  activity,  etc.  A  current  proposed  model  assumes  a 
constant  working  temperature  [1],  which  does  not  incorpo¬ 
rate  the  reality  that  this  working  temperature  dynamically 
varies  at  different  load  currents,  as  well  as  during  startup  and 
shut-down  of  the  fuel  cell  system.  Some  authors  proposed 
improved  models,  with  Amphlett  et  al.  [2]  using  the  first  em¬ 
pirical  thermal  model,  and  Gurski  et  al.  [3]  considering  the 
reactant  flows  and  coolant  control  based  upon  the  previous 
model.  Others  proposed  models  calculating  the  temperature 
variation  of  the  stack,  cell  [4-10]  or  two  electrodes  and  ME  A 
[11,12].  B.  Wetton  et  al.  [13]  proposed  an  explicit  thermal 
model  to  analyze  the  temperature  gradient  of  different  layers 
in  the  fuel  cell  stack  considering  the  stack  asymmetric  effects, 
which  does  not  include  dynamics.  Recently,  M.  Sundaresan 
published  the  most  detailed  ID  thermal  dynamic  model  [14]. 
However,  the  flow  of  species  at  the  inlet  must  be  the  same 
as  that  at  the  outlet.  Thus,  no  fluid  dynamics  is  considered  in 
the  model. 

Secondly,  the  proton  transport  in  the  membrane  and  its 
associated  ohmic  losses  mainly  determine  the  characteristics 
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Nomenclature 

Alphabets 

a 

species  activity 

A 

area  (m2) 

C 

mass  concentration  (kg  m-3) 

D 

diffusion  coefficient  (m2  s-1) 

F 

faraday  number 

G 

gibbs  free  energy  (J  mol-1) 

H 

enthalpy  (J  mol-1) 

i 

current  density  (Am-2) 

j 

exchange  current  density  (Am-2) 

l 

thickness  (m) 

m 

mass  (kg) 

M 

mole  mass  (kg  mol-1) 

nd 

electro-osmotic  drag  coefficient 

N 

mole  flux  (mol  s-1  m-2) 

P 

pressure  (partial  pressure)  (Pa) 

R 

universal  gas  constant 

Rmem 

proton  transfer  resistance  (£2) 

^ab 

electrical  resistance  (Q) 

s 

entropy  (J  mol-1  K-1) 

T 

temperature  (K) 

W 

mass  flux  (kg  s-1  m-2) 

Greek  symbols 

£ 

porosity 

water  uptake  coefficient 

p 

density  (kg  m-3) 

X 

tortuosity 

Superscripts  and  subscripts 

an 

anode 

ca 

cathode 

cv 

control  volume 

d 

gas  diffusion  layer 

g 

gas 

i 

index 

1 

liquid 

mem 

membrane  layer 

ref 

reference  value 

sat 

saturation 

sou 

source 

of  ohmic  polarization.  The  proton  conductivity  has  been  re¬ 
garded  as  constant,  temperature  dependent  [1]  or  temperature 
and  water  concentration  dependent  variables  [15].  Recently, 
Pukrushpan  et  al.  [16]  proposed  the  most  comprehensive 
model  that  considers  the  dependence  of  the  proton  conduc¬ 
tivity  on  the  water  concentration  and  temperature.  However, 
the  water  concentration  of  the  membrane  is  obtained  from  the 
membrane  relative  humidity  (RH)  on  an  average  of  the  an¬ 
ode  and  cathode  RH.  In  fact,  the  RH  in  the  anode  and  cathode 


varies  rapidly,  while  the  RH  in  the  membrane  does  slowly  be¬ 
cause  the  amount  of  water  residing  in  both  sides  is  relatively 
less  than  in  the  membrane  [15]. 

Thirdly,  the  oxygen  concentration  in  the  GDL  on  the  cath¬ 
ode  side  is  continuously  changing  in  operating  environments 
and  significantly  affects  the  performance  of  the  cell.  There¬ 
fore,  plenty  of  models  considering  multi-phase  multi-species 
have  been  employed  to  investigate  the  transport  phenom¬ 
ena  in  the  GDL.  However,  those  models  do  not  consider  the 
dynamics.  Recently,  Pukrushpan  et  al.  proposed  a  dynamic 
model  with  lumped  parameters  to  predict  the  gas  dynamics 
in  a  cathode  electrode,  which  does  not  consider  the  effects  in 
the  GDL  [16].  In  this  paper,  we  use  a  ID  single-phase  model 
to  represent  the  dynamics  present  in  the  GDL. 


2.  Model  setup  and  assumptions 

The  model  has  been  developed  on  the  basis  of  layers  in  a 
cell  that  consist  of  a  MEA,  two  gas  diffusion  layers  and  two 
gas  channels  sandwiched  by  two  coolant  channels,  as  shown 
in  Fig.  1.  The  input  variables  for  the  model  are  current  load, 
mass  flow  rate,  the  gas  components  fraction,  temperature, 
pressure  and  relative  humidity  of  reactants  as  well  as  the 
temperature  and  velocity  of  coolants  at  the  inlets. 

The  main  assumptions  made  for  the  new  model  are  as 
follows: 

1.  Reactants  are  ideal  gases. 

2.  There  is  no  pressure  gradient  between  the  anode  and  cath¬ 
ode  side;  it  means  no  convection  but  only  diffusion  for  gas 
transport  is  considered. 

3.  There  is  no  gas  pressure  drop  from  the  inlet  to  the  outlet 
of  the  gas  channel. 

4.  The  temperature  gradient  is  linear  across  the  layers  in  a 
fuel  cell. 

5.  The  thermal  conductivity  for  the  materials  in  a  fuel  cell  is 
constant. 

6.  There  is  no  contact  resistance. 

7.  Anodic  over-potential  is  negligible. 


Fig.  1.  Schematic  simulation  domain. 
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8.  There  is  no  current  density  gradient  across  the  cathode 
catalyst  layer;  it  means  that  the  reactants  completely  re¬ 
acted  as  soon  as  it  reaches  to  the  cathode  catalyst  layer 
surface. 

9.  Based  on  these  assumptions,  five  sub-models  have  been 
developed  and  are  described  in  the  following  sections. 


3.  Model  description 

3.1.  Electrochemical  model 


changes  in  a  controlled  volume  equals  the  sum  of  the  energy 
exchange  at  boundaries  and  internal  energy  resources.  In  fact, 
the  energy  exchanges  at  boundaries  occur  by  three  factors: 
(a)  the  mass  flow  into  each  volume;  (b)  the  conduction  heat 
transfer  across  the  cell;  (c)  the  convection  heat  transfer  occur¬ 
ring  between  bipolar  plates  with  the  coolant  and  the  reactants. 
Thus,  the  thermal  dynamic  behavior  can  be  described  with 
the  following  energy  conservation  equation: 


]T  cPic 


z,mass 


drcv 

d  t 


i 


Generally,  the  overall  chemical  reaction  of  the  PEM  fuel 
cell  can  be  described  by  using  the  following  expressions, 
illustrating  that  a  chemical  reaction  of  hydrogen  and  oxy¬ 
gen  molecules  produces  electricity,  water  and  heat  as  a  by¬ 
product: 

H2  +  \02  H20  +  Q  res  T  Tceii 

The  output  cell  voltage  Vceii  is  the  difference  between  the 
open  circuit  voltage  (OCV)  Eq  and  over-potentials  77  and 
Vohm- 


^  ^  ril'mAcQ\\Cp j  (7]n 

Tcv)  T 

GConvAcell 

s 

Vs - X/ - ^ 

mass  flow  in 

convection  heat  transfer 

T  0CondACell 
^  „  ✓ 

T  j3sou 

conduction  heat  transfer 

sources 

On  the  other  hand,  the  internal  energy  source  is  composed  of 
the  entropy  loss  and  the  chemical  energy  required  for  protons 
to  overcome  the  barrier  of  the  over-potentials  in  both  catalyst 
layers  (Eq.  (7)).  In  addition,  other  heat  sources  are  ohmic 
losses  caused  by  a  transport  of  electrons  and  protons  in  the 
cell: 


By  neglecting  the  dependence  of  the  OCV  on  the  reactant 
pressure,  the  relationship  between  the  OCV  and  the  tempera¬ 
ture  can  be  simplified  with  the  empirical  parameter  dEo/dT.  If 
the  reactant  is  ideal,  its  activity  can  be  described  by  using  Eq. 
(2),  where  index  i  indicates  H2  and  O2,  while  P(  is  the  partial 
pressure  of  gas  components,  and  Po  is  the  overall  pressure  of 
both  the  anode  and  cathode  side.  Then,  Eo  can  be  derived  by 
modifying  the  Nernst  Eq.  (3). 


A) 


Q  sou  —  ^cell  ( 


TAS 


+  0  +  iAcqwRqIq  J 


In  fact,  the  change  of  entropy  due  to  the  electrochemical  re¬ 
action  (Eqs.  (8)  and  (9))  in  both  of  the  catalysts  sides  pre¬ 
dominantly  influences  the  energy  sources  term  according  to 
the  calculation  shown  below. 


«2  ^  2H+q)  +  2e(pt)-  (8) 

02  +  2H+  +  2e~  =  H20(i)  (9) 


Eo  =  £ref  +  ~^(T  -  7W)  +  —  ln(<4<5)  (3) 

The  anodic  over-potential  is  negligible;  while  the  ij  repre¬ 
sents  the  over-potential  of  the  cathode  catalyst  layer.  Under 
the  further  assumptions  that  the  asymmetric  parameter  of  the 
reaction  is  (1)  and  (8),  the  Butler- Volmer  equation  leads  to 
Eq.  (4)  that  describes  the  over-potential  on  the  cathode  side. 


i  =  Jo 


Acata,eff 
A  cell 


Po2  [H+]  1 

P02,ref  [H+]ref  ^  \RTJ 


In  order  to  obtain  the  entropy  change  of  these  reactions,  the 
zero  point  of  semi-absolute  entropy  is  taken  as  a  reference 
according  to  [17]: 

•v|H+q)]  =  0  (10) 

The  entropy  of  an  electron  obtained  from  the  standard  hydro¬ 
gen  electrode  results  in  the  following  equations  [17]: 


ASshe  = 


AT/she  -  AGshe 
T 


(11) 


The  ohmic  over-potential  Vohm  is  determined  by  the  product 
of  the  current  density  and  the  proton  resistance  7?mem  accord¬ 
ing  to  Ohm’s  law  (5). 

Eohm  —  IP  mem  (5) 


1  1  1 

j[e(pt)“]  =  -s[  H2]  =  65.29  JmoU1  K-1  (12) 

2 

Therefore,  the  entropy  change  of  the  cathode  reaction  is  equal 
to  the  sum  of  that  in  water,  oxygen  and  electron: 


3.2.  Thermal  model  for  a  cell 

If  a  cell  is  assembled  with  cubic  layers  whose  thermo¬ 
physical  properties  are  isotropic  and  constant,  and  then  ac¬ 
cording  to  the  energy  conservation  equation,  the  total  energy 


ASCa  =  s[H20(i)]  -  ^[02]  -  2s[e(pt)  ]  =  69.91 


1 

2 


X  205.03  -  2  X  65.29  =  -163  J  mol-1  K"1 


(13) 
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If  the  anode  is  assumed  as  a  standard  electrode,  the  anodic 
entropy  change  becomes  0. 


3.3.  Proton  conducting  model  for  membrane 


The  membrane  resistance  is  a  function  of  the  temperature 
and  water  content  in  a  membrane  layer,  which  is  described 
as  follows  [16]: 


^mem  — 


■  mem 


{bn  ^mem  £12) exp  (b2  (3^3 


mem 


(14) 


where  the  temperature  Tmcm  can  be  derived  from  the  previ¬ 
ous  Eq.  (6),  while  the  membrane  water  content  Amem  can  be 
described  by  using  the  water  mass  concentration  [15]  and  the 
mass  conservation  equation  [16]: 


A 


mem 


C  h2o  ,mass  /Mh2o 
Airy, mem/ Mmcm  ^CH20,mass/ 


(15) 


dewater  d(Cn20,  mass cell  ^  mem) 


d  t 


d  t 

Wele,mem,an  -  Wele  ,mem,ca  + 
VEdiff  ,mem,ca  +  Wdiff  ,  mem,  an 


(16) 


The  electro-osmotic  driving  force  created  by  the  different 
electrochemical  potential  at  the  anode  and  cathode  deter¬ 
mines  the  water  mass  flows  of  Wele,mem,an  and  Wele,mem,ca 
at  the  boundaries  of  the  membrane  layer.  In  addition,  the 
diffusion  caused  by  the  water  concentration  gradient  at  the 
two  boundaries  makes  up  the  mass  flows  of  Wdiff, mem, an 
and  Wdiff,mem,ca-  Those  relationships  are  described  by  Eqs. 
(17)— (19),  proposed  by  Spriger  [24]. 


n  d 


=  0.0029A^em  +  0.05/,mem  -  3.4  x  1(T19 


^ele.mem,;  —  ^Vvater -Veiled,! 

r 


Wdiff  ,mem,i  —  ^water^.  cell  ^  water 


Ci  Emid 


l 


mem 


(17) 

(18) 

(19) 


Hence,  the  diffusion  coefficient  Dwa ter  and  the  water  concen¬ 
tration  Q  are  calculated  from  the  empirical  Eq.  [24]  expressed 
as  a  function  of  membrane  water  content  Xmem: 


D- water  —  7)(Amem)  exp 


1 


Tmem 


(20) 


77  (7- mem)  —  ^ 


10"6 

2  >  7-mem 

10-6(1  +  2(Amem  —  3)) 

3  2  7.mcm  >  2 

10-6(3  —  1.67(7.mem  —  3)) 

4.5  >  Amem  >  3 

1.25  x  10“6 

7-mem  :2  4.5 

(21) 

The  boundary  water  content  A/  is  a  function  of  water  activity 
au  which  is  calculated  from  the  water  vapor  partial  pressure: 


'  0.043  +  17.81a,-  -  39.85 a]  +  36 a] 
14+1.4 (at  -  1) 

16.8 


1  >  a[  >  0 
3  >  ci  [  >  1 
3  <  a[ 

(22) 


(23) 


3.4.  Proton  dynamic  model  in  the  cathode  catalyst  layer 


The  dynamic  behavior  of  a  fuel  cell  at  a  load  is  investigated 
by  experiments.  When  the  output  current  changes  abruptly, 
the  output  voltage  of  the  fuel  cell  reacts  with  an  overshoot 
[18].  These  dynamics  result  from  different  physical  phenom¬ 
ena  of  reactants  and  their  chemical  reaction  in  the  cell,  such 
as  dynamics  filling  in  the  gas  flow  channel,  diffusing  reac¬ 
tants  through  the  GDL  and  reacting  process  in  the  double 
layer  at  the  interface  of  electrodes.  Ceraolo  et  al.  explained 
the  dynamic  effects  with  a  relationship  between  the  number 
of  mobile  protons  and  water  content  [1].  As  a  matter  of  fact, 
when  the  current  density  increases,  the  hydration  of  the  poly¬ 
meric  electrolyte  near  the  cathode  catalyst  tends  to  rise  as 
well;  consequently,  the  proton  concentration  near  the  cath¬ 
ode  catalyst  increases  rapidly.  On  the  other  hand,  the  proton 
concentration  will  decrease  slowly  at  a  decrease  of  current. 
Thus,  the  dynamics  can  be  described  by  the  following  differ¬ 
ential  equation  using  the  proton  concentration  as  a  variable 
[1]: 


1  +  aH+/3 
th+ 


(24) 


CH+  =  [H+]/[H+  ]ref  is  the  dimensionless  proton  concentra¬ 
tion,  <$( )  the  Heaviside  function,  and  rH+  and  aH+  are  empir¬ 
ical  parameters. 

Fig.  2  shows  the  calculated  response.  The  voltage  de¬ 
creases  quickly  when  the  current  density  increases.  However, 
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Fig.  2.  Voltage  response  by  a  consideration  of  proton  concentration. 
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the  voltage  first  reaches  its  highest  value  and  then  damps  with 
a  time  constant  that  is  associated  with  the  proton  concentra¬ 
tion,  as  the  current  density  decreases. 

3.5.  GDL  reactant  model  for  cathode 


■  i 


.  Current  Density! . * 

5  6 

Time  (sec) 


Air  contains  not  only  oxygen  but  also  nitrogen  and  wa¬ 
ter  vapor.  The  air  entering  the  cell  diffuses  through  the 
GDL  before  reaching  the  catalyst  layer.  The  diffusion  ef¬ 
fect  is  described  by  using  the  mass  conservation  (25)  and 
Stefan-Maxwell  equations  (26): 


£g  dP[  dNi 
~RT~df  +  ~dy 


£  g  3  Pi 

x1  3  y 


£  - PM> 


(25) 

(26) 


Hence,  i,  k  e  (1,  3),  where  P\  is  the  oxygen  partial  pressure, 
and  P2  =  T>sat(7)  and  P3  are  the  water  vapor  and  the  nitro¬ 
gen  partial  pressure,  respectively.  The  diffusion  coefficient 
P caPik  =  A*,e ff  =  Afc,eff and  the  cathode  pressure  of  FCa  is 
the  summation  of  the  species  partial  pressures.  The  parameter 
r  is  a  constant  describing  the  pore  curvature  of  the  GDL. 

The  partial  differential  equation  (PDE)  systems  above  can 
further  be  simplified  by  using  the  following  PDE  [  1  ] ,  whereby 
§  is  the  dimensionless  distance  y/l&: 


dPo2  d2  Po2  i  dPo2 

- -  =  CO - ^  f - - 

3 1  df  4  F  3§ 

1 

T2l2d((Psat/Dn)  +  (Pea  -  P^)/Du) 


(27) 

(28) 


RT 

£g^d(^ca  —  ^sat) 


(29) 


In  order  to  investigate  the  effects  of  the  GDL  on  dynamics, 
simulations  have  been  run  to  analyze  the  relationship  between 
the  GDL  thickness  and  the  dynamics  of  oxygen  by  diffusions. 
Fig.  3  shows  the  dynamic  response  of  the  oxygen  partial  pres¬ 
sure  at  the  interface  of  the  cathode  GDL  and  cathode  cata¬ 
lyst  layer.  The  results  show  that  the  oxygen  partial  pressure 
drops  rapidly  when  a  step  current  (Fig.  4)  is  applied.  Thus, 


Fig.  4.  Current  input. 


the  concentration  over-potential  increases.  Accordingly,  the 
thickness  of  the  GDL  is  one  of  the  factors  influencing  the 
dynamic  response.  Moreover,  the  steady  state  is  reached  by 
the  thin  GDL  more  quickly  than  by  the  thick  one. 

Further  analysis  has  been  undertaken  to  discover  how  the 
reactant  partial  pressure  is  distributed  along  the  GDL  by  the 
given  pressure  (Fig.  5)  and  current  step.  Fig.  6  shows  the  cath¬ 
ode  oxygen  partial  pressures  and  their  responses  depending 
upon  the  thickness  ratio.  The  analysis  shows  that  the  dynamic 
response  of  the  oxygen  partial  pressure  is  highly  dependent 
upon  the  geometrical  locations.  When  the  cathode  inlet  pres¬ 
sure  changes,  the  pressure  at  the  catalysts  side  responds  with 
a  time  delay  before  it  has  reached  the  steady  state,  which 
is  caused  by  the  diffusion  of  the  reactant.  Accordingly,  the 
over-potential  cannot  be  manipulated  instantly. 

The  dynamic  responses  of  the  oxygen  concentration  at  the 
catalyst  layer  are  illustrated  in  Fig.  7.  The  oxygen  concen¬ 
tration  is  strongly  influenced  by  the  thickness  of  the  GDL. 
The  thinner  the  layer  becomes,  the  shorter  the  response  time 
gets.  On  the  other  hand,  when  the  inlet  pressure  increases 
(Fig.  5),  the  partial  pressure  at  the  catalysts  tends  to  follow 
its  increase,  but  the  amounts  of  the  recovered  partial  pressure 
compared  to  Fig.  3  depend  on  the  thickness.  Therefore,  the 
settle  times  to  the  steady  state  become  constant  regardless  of 
the  thickness  of  the  GDL. 


Fig.  5.  Dynamic  cathode  pressure  input. 
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Fig.  3.  Oxygen  concentration  response  dependent  upon  GDL  thickness  at 
P=  1  bar  and  T  =  353  °K. 
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Fig.  6.  Oxygen  concentration  response  to  the  dynamic  current  input  and 
cathode  pressure  at  different  depths. 
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Fig.  7.  Oxygen  concentration  dynamic  response  for  different  GDL  thick¬ 
nesses. 

In  this  part,  the  effects  of  material  properties  on  the  dy¬ 
namics  are  analyzed  by  a  given  current  step.  Figs.  8  and  9 
show  the  dependences  on  porosity  and  tortuosity.  Generally, 
the  GDL  with  a  high  porosity  allows  less  pressure  drops  than 
that  with  the  lower  one.  However,  the  GDL  with  a  high  tortu¬ 
osity  causes  higher  pressure  drops  than  the  low  one  because 
of  a  long  path  for  the  oxygen  transported.  When  a  step  current 
is  applied  to  the  cell,  the  oxygen  consumption  on  the  catalysts 
side  will  be  increasing.  Instantly,  the  high  porosity  enables 
more  oxygen  to  transfer  from  the  inlet  to  the  catalysts  side, 
and,  subsequently,  the  oxygen  partial  pressure  at  the  catalysts 
can  quickly  follows  the  pressure  changes  at  the  inlet.  Thus, 
the  GDL  with  a  high  porosity  dynamically  responds  to  the 
pressure  increase.  When  the  tortuosity  increases,  the  dynamic 
response  time  slows  from  the  same  effects  as  prolonged  ge¬ 
ometry  and  the  associated  pressure  drop.  The  settle  times 
remain  unchanged,  as  in  the  analysis  for  different  thickness. 


Time  (sec) 

Fig.  8.  Dynamic  response  of  oxygen  concentration  for  different  porosities. 


Fig.  9.  Dynamic  response  of  oxygen  concentration  for  different  tortuosities. 


4.  Simulation  results 

All  of  the  aforementioned  models  are  coded  with  MAT- 
LAB  and  C.  Multi-run  simulations  have  been  conducted  to  in¬ 
vestigate  the  static  and  dynamic  behavior  of  a  single  cell.  The 
static  behavior  is  analyzed  by  calculating  the  typical  polar¬ 
ization,  which  takes  into  account  variables  such  as,  working 
temperature,  pressure,  stoichimetric  number  and  relative  hu¬ 
midity  (RH)  at  the  cathode  inlet.  The  dynamic  characteristic 
considers  two  aspects,  the  startup  and  the  transient  response 
on  the  current  as  a  step  load. 

4.1.  Parameters  and  reference  data 

The  parameters  and  reference  data  for  the  models  chosen 
are  as  follows  (see  Table  1),  and  they  are  partially  empirical 
[1,16,21]. 

4.2.  Static  behavior 

Fig.  10  shows  the  temperature  dependent  I-V  characteris¬ 
tics  from  333  to  363  K  with  a  step  of  10  K.  As  the  temperature 
increases,  the  water  removal  will  be  more  eased.  The  effects 
are  considerably  high  at  the  range  of  the  higher  cell  current, 
where  more  water  is  produced.  This  result  is  comparable  with 
the  CFD  analysis  [19]. 

Fig.  1 1  shows  the  pressure  dependent  polarization  curve. 
As  the  pressure  increases,  the  oxygen  concentration  at  the 
surface  of  the  catalysts  layer  tends  to  increase,  too.  Thus,  the 
concentration  over-potential  gets  lower.  Otherwise,  the  over¬ 
potential  becomes  higher  because  of  the  oxygen  starvation. 

Table  1 


Empirical/reference  parameters 


Electronic  reactions  model 

A) 

1 .0  bar 

Lcf 

343.15  K 

Act 

0.975  V  [1] 

d  E0/dT 

0.00027  V  K-1  [1] 

^cata.cffA  cell 

/(/,  T,  P02)  [1] 

Gas  transport  model 

f(P,  7)m2  s-1  [1] 

P  sat 

f(T)  bar  [1] 

Thermal  model 

Hgas 

f(P,  T)  [21] 

Cp-gas 

f(P,  T)  [21] 

P  gas 

f(P,  T)  [21] 

Proton  conducting  model 

b  n 

0.5139  [16] 

b\2 

0.326  [16] 

hi 

350  [16] 

h 

0.0126  [16] 

nd 

/(CWater)  [16] 

£>w 

f(T,  CWater)  [16] 

Proton  concentration  model 

aH+ 

5.87E-12(m2  A-1)3  [1] 

th+ 

12.78  s  [1] 
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Fig.  10.  I-V  curve  for  different  cell  working  temperature.  Cell:  P=  1.0  bar. 


Fig.  12.  I-V  curve  for  different  stoichimetric  number.  Cell:  T  =  353. 15  K, 
P  =  1 .0  bar. 


Fig.  11.  I-V  curve  for  different  cell  gas  pressure.  Cell:  T=353.15K. 

The  results  are  also  comparable  with  the  CFD  analysis 
[19]. 

Fig.  12  shows  the  stoichimetric  number  dependent  I-V 
curve  at  a  constant  temperature  and  pressure.  When  the  sto¬ 
ichimetric  number  is  low,  the  removal  of  water  at  the  cathode 
outlet  flow  decreases.  Thus,  the  water  concentration  in  the 
membrane  layer  increases.  Consequently,  the  membrane  re¬ 
sistance  and  the  resulting  ohmic  over-potential  become  lower. 
However,  the  low  stoichimetric  number  adversely  affects  the 
cathode  over-potential  at  the  high  current  because  of  the  ex¬ 
cessive  water  in  the  catalysts. 

Fig.  13  shows  the  output  voltage  of  a  cell  influenced  by 
relative  humidity  at  the  cathode  inlet.  When  the  humidity  in¬ 
creases  at  the  cathode  side,  the  air  transported  to  the  catalysts 


Fig.  13.  I-V  curve  for  different  RH  at  the  cathode  inlet.  Cell:  T  =  353.15  K, 
P  =  1 .0  bar. 

will  be  blocked,  and  the  cathode  over-potential  will  increase, 
especially  at  the  high  current  range.  The  results  are  compa¬ 
rable  to  the  experimental  data  [20] . 

4.3.  Dynamic  behaviors 

For  a  startup  of  dynamic  simulations,  initial  values  are 
necessary  for  variables  such  as  layer  temperature,  membrane 
water  concentration,  GDL  air  and  oxygen  concentration  and 
gas  channel  pressure.  Fig.  14  shows  the  setup  for  the  simula¬ 
tion  of  a  single  cell  from  layer  1  to  11. 

Geometrical  and  thermo-physical  data  for  the  layers  are 
summarized  in  Table  2. 


Table  2 

Simulation  data 


Thickness  (m) 

Heat  conductivity  (Wm  1  k  1 ) 

Heat  capacity  (J  Kg  1  K  1 ) 

Density  (Kg  m  3 ) 

GDL 

0.0004 

65 

840 

2000 

Catalyst  layer 

0.000065 

0.2 

770 

387 

Membrane  layer 

0.000183 

0.21 

1100 

1967 

Gas  channel 

0.001 

52 

935 

1400 

Plate 

0.001 

52 

935 

1400 

Coolant  channel 

0.001 

30 

935 

1400 

GDL  porosity 

0.5 

GDL  tortuosity 

3.725 

Bipolar  plate  contact  area  percentage 

0.55 

Membrane  molecular  mass 

1 . 1  Kg  mol- 1 

Fuel  cell  area 

0.0367  m2 

Fuel  cell  active  area 

0.03  m2 
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Cathode  Bipolar  Plate 


Anode  Bipolar  Plate 


Fig.  14.  Simulation  setup. 


4.3.1.  Startup 

The  startup  temperature  for  the  cell  model  is  initially  set 
to  298.15  K.  The  value  of  current  density  is  increased  con¬ 
tinuously  for  the  first  350  s  in  order  to  quickly  raise  the  tem¬ 
perature  to  353.15  K,  which  is  assumed  as  a  typical  work¬ 
ing  temperature.  In  addition,  a  temperature  controller  is  built 
in  the  simulation,  as  if  a  coolant  subsystem  turned  on  and 
off  at  this  set  point  to  extract  the  excessive  heat  produced 
in  the  cell.  Fig.  15  shows  the  dynamic  behavior  of  the  tem¬ 
perature  for  different  layers  and  voltage,  as  well  as  the  ef¬ 
ficiency  during  a  startup.  It  took  8  min  for  the  cell  to  reach 
to  the  working  temperature  (Fig.  15a).  Generally,  the  tem¬ 
perature  profiles  in  each  of  the  layers  tend  to  follow  the 
current  waveform,  because  of  the  associated  energy  losses 
occurring  in  the  layers.  Particularly,  the  temperature  in  the 
membrane  and  catalysts  layer  is  the  highest,  which  results 
from  ohmic  losses  due  to  the  membrane  resistance  and  the 
heat  released  by  the  chemical  reaction.  The  average  differ¬ 
ence  of  temperature  between  these  two  layers  and  other  lay¬ 
ers  on  the  anode  and  cathode  side  amounts  to  3  and  2K, 
respectively.  Corresponding  voltage  and  power  are  shown 
in  Fig.  15c.  When  the  current  increases,  the  over-potentials 
increase,  and,  subsequently,  the  voltage  and  power  decrease. 
While  the  temperature  is  rising,  the  voltage  fluctuates  slightly 
during  the  startup.  When  the  current  had  been  kept  constant 
after  the  350  s,  the  amount  of  water  generated  in  the  cat¬ 
alyst  layer  becomes  constant.  On  the  other  hand,  the  con¬ 
tinuously  increased  temperature  leads  to  a  high  saturation 
pressure  in  the  cell,  which  enables  water  residing  in  the  cat¬ 
alysts  to  be  quickly  removed  [23].  Otherwise,  water  would 
be  flooding  and  blocking  further  influx  of  the  oxygen  into 
the  catalysts.  Therefore,  the  cell  voltage  increases  rapidly. 
Thereafter,  the  water  concentration  in  the  membrane  contin¬ 
uously  decreases  by  the  electro-osmotic  force  and  diffusion 
effects,  and  the  corresponding  proton  conductivity  will  be 
decreased.  Thus,  the  cell  voltage  slightly  drops  after  the  tem¬ 
perature  has  reached  a  steady  state.  The  overall  efficiency  of  a 
cell  is  also  strongly  influenced  by  the  variation  of  temperature 
(Fig.  15d). 

Fig.  16  shows  the  dynamic  behavior  of  the  temperature 
distribution  across  the  fuel  cell  at  50  s  and  7  min  during  a 
startup  process.  The  catalyst  on  the  cathode  side  shows  the 
highest  peak  because  of  the  losses  associated  with  the  over¬ 
potential  being  higher  than  on  the  anode  side.  For  the  first 
30  s,  the  temperature  rises  slowly  because  of  the  slow  slope  of 
the  current.  The  maximum  difference  of  temperature  between 
the  GDL  at  the  anode  side  and  the  membrane  has  been  shown 
to  reach  1  K.  The  higher  the  current  is,  the  more  heat  is  gener- 


Dynamic  temperature  response  of  different  layer  of  single  cell 


Single  cell  efficiency 
Fig.  15.  Simulation  results. 

ated.  Thus,  the  rise  of  temperature  in  the  layers  is  accelerated. 
When  the  startup  is  ended  after  7  min,  the  temperature  in  the 
catalyst  rises  up  to  353  K,  approximately  2  K  higher  than  the 
anode  side  catalyst.  Then,  the  peak  point  of  the  temperature 
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Fig.  16.  Temperature  gradient  across  the  cell  (left  to  right:  cathode  coolant 
channel  to  anode  coolant  channel):  (a)  50  s  and  (b)  7  min  after  startup. 


is  moved  from  the  catalyst  to  the  membrane,  which  results 
from  the  dehydration  of  the  membrane  and  the  associated 
increase  of  losses.  The  dehydration  is  mainly  caused  by 
diffusions  of  water  from  the  membrane  to  both  sides  because 
of  higher  water  concentration  in  the  membrane  than  in  the 
gas  channel  sides.  On  the  other  hand,  the  increased  number  of 
protons  transported  takes  up  more  water  from  the  membrane 
to  the  cathode.  Consequently,  the  resistance  in  the  membrane 
is  increased  and  shows  the  highest  temperature  among  the 
layers. 


4.3.2.  Transient  response 

In  order  to  analyze  the  dynamic  response  on  a  power  de¬ 
mand,  a  step  current  with  0. 8-0.4  A  cm-2  at  the  600  s  is  ap¬ 
plied.  Fig.  17  shows  the  response  of  the  temperature  in  the 
different  layers  (b),  voltage  (c),  power  output  (c)  and  effi¬ 
ciency  (d). 

The  operating  temperature  is  automatically  controlled  by 
the  coolant  system,  the  reference  value  for  which  is  set  to 
353.15  K.  The  on-off  control  of  the  coolants  causes  a  slight 
fluctuation  of  the  temperature  waveforms  until  600  s.  When 
the  current  suddenly  decreases,  the  heat  generated  at  the  cath¬ 
ode  catalyst  and  membrane  layer  decreases  rapidly  and  leads 
to  a  temperature  drop  at  these  two  layers.  Then,  the  coolant 
system  is  turned  off.  The  heat  is  transferred  by  the  tempera¬ 
ture  gradient  from  the  layers  into  the  bipolar  plates  and  stored 
there.  Thus,  the  temperature  of  both  bipolar  plates  tends  to 
increase,  and  the  temperature  gradient  begins  to  decrease.  As 
a  result,  the  amount  of  heat  removed  from  the  catalyst  and 
membrane  layers  is  again  decreased. 
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Fig.  17.  Analysis  of  transient  behavior  of  temperature,  voltage,  power  and 
efficiency  upon  a  current  step. 


When  the  temperature  of  the  catalyst  and  the  membrane 
layer  reaches  its  lowest  point,  the  temperature  of  all  the  layers 
will  rise  again  due  to  the  accumulated  heat  after  the  coolant 
is  turned  off.  Finally,  it  reaches  the  steady  state  after  around 
10  s. 

The  effects  of  the  temperature  variation  on  the  output 
voltage  are  slightly  different  from  the  temperature  profiles. 
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When  the  current  decreases,  the  voltage  rapidly  increases 
and  shows  a  slight  overshoot.  Then,  the  voltage  slowly  el¬ 
evates  and  drops  back  to  a  steady  state.  The  first  overshoot 
of  the  voltage  results  from  the  variation  of  the  proton  con¬ 
centration  in  the  cathode  catalyst  layer,  while  the  first  half 
of  the  voltage  arc  is  caused  by  the  temperature  increase  in 
all  the  layers,  and  the  rest  by  the  decrease  of  water  concen¬ 
tration  in  the  membrane  layer.  The  efficiency  profile  is  the 
same  as  the  behaviors  of  the  output  voltage,  but  the  power 
remains  almost  the  same  as  before  the  current  was  applied. 
The  voltage  response  is  comparable  to  an  experimental  result 
[22]. 


5.  Conclusion 

In  this  paper,  a  new  dynamic  model  for  the  PEM  fuel 
cell  is  proposed  to  predict  responses  of  the  electric  load  on 
the  cell.  Emphasis  is  placed  on  the  temperature  response 
to  the  dynamic  load.  The  model  is  constructed  with  layers 
that  include  membrane,  catalysts,  gas  diffusion  layers  and 
bipolar  plates.  The  models  for  the  layers  are  separately  and 
mathematically  described.  The  analyses  delivered  surpris¬ 
ing  results  in  the  effects  of  electric  loads  on  temperature, 
voltage  and  efficiency.  Particularly,  description  of  interre¬ 
lated  physical  variables  by  varying  temperature  provides  a 
more  realistic  tool  that  can  be  utilized  for  deriving  design 
parameters  of  a  fuel  cell  stack  and  systems.  As  an  example, 
a  startup  process  and  the  associated  effects  are  analyzed  in 
detail. 


6.  Future  work 

In  the  near  future,  the  model  will  be  integrated  into 
the  exiting  system  model  and  implemented  on  a  RT  sys¬ 
tem,  which  enables  us  to  optimize  a  fuel  cell  powered 
system  and  design  advanced  controls.  Further  improve¬ 
ments  are  anticipated  by  taking  into  account  the  gas 
dynamics  in  the  flow  channel.  At  the  same  time,  we 
plan  to  validate  the  models  in  a  close  collaboration  with 
industry. 
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